home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Monster Media 1996 #15
/
Monster Media Number 15 (Monster Media)(July 1996).ISO
/
prog_c
/
cuj0696.zip
/
DWYER.ZIP
/
SPECTRAL.TST
/
SPECTRAL.CX
< prev
next >
Wrap
Text File
|
1996-03-23
|
8KB
|
318 lines
/* ============ */
/* spectral.cx */
/* ============ */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <xtendefs.h>
#include <specdefs.h>
#include <calcmert.cx>
#include <euclstep.cx>
#include <getctlda.cx>
#include <getgcd.cx>
#include <prtspdat.cx>
#include <putxflt.cx>
#include <putxfrac.cx>
#include <putxint.cx>
#include <xdot.cx>
#include <xmatprt.cx>
/* ==================================================================== */
/* MAIN - Performs Spectral Test on multiplier of CRNG */
/* ==================================================================== */
void
main()
{
short HowMany, InputOK, MaxDimen;
USHORT Modulus[NE], Mult[NE * 2];
USHORT Tmp1[NE], Tmp2[NE];
InputOK = GetControlData(Mult, Modulus, &MaxDimen, &HowMany);
if (InputOK)
{
USHORT XH[NE], XHP[NE], XP[NE], XPP[NE],
XR[NE], XS[NE], XU[NE], XV[NE];
USHORT MSqd[NE]; /* Modulus Squared */
USHORT XMerit[NE]; /* Holds Accuracy Results */
printf("\n%22d%41s", HowMany, "Number of Multipliers");
printf("\n%20s", "");
PutXInt(Mult, 20);
printf(" Multiplier");
if (HowMany == 1)
{
/* Do Nothing */
}
else
{
/* ------------------------------- */
/* Print Label on First Multiplier */
/* ------------------------------- */
printf(" # 1\n%20s", "");
/* ----------------------- */
/* Print Second Multiplier */
/* ----------------------- */
PutXInt(Mult + NE, 20);
/* ----------- */
/* Print Label */
/* ----------- */
printf(" Multiplier # 2");
}
printf("\n%20s", "");
PutXInt(Modulus, 20);
printf(" Modulus\n%21s", "");
printf("%-19d Maximum Dimension of Test\n", MaxDimen);
if (HowMany == 1)
{
/* -------------- */
/* S1: Initialize */
/* -------------- */
XCOPY(Mult, XH); /* h <- a */
XCOPY(Modulus, XHP); /* m <- h' */
XCOPY(EONE, XP); /* p <- 1 */
XCOPY(EZERO, XPP); /* p' <- 0 */
XCOPY(Mult, XR); /* r <- a */
XSQR(Mult, XS);
XADD(XS, EONE, XS); /* S <- 1 + a^2 */
/* ------------------ */
/* S2: Euclidean Step */
/* ------------------ */
EuclStep(XH, XHP, XP, XPP, XS, XU, XV);
/* -------------- */
/* S3: Compute V2 */
/* -------------- */
XSUB(XU, XH, XU); /* u -= h */
XSUB(XV, XP, XV); /* v -= p */
XSQR(XU, Tmp1); /* u^2 */
XSQR(XV, Tmp2); /* v^2 */
XADD(Tmp1, Tmp2, Tmp2); /* u^2 + v^2 */
if (XGT(XS, Tmp2)) /* if(s > Tmp2) */
{
XCOPY(Tmp2, XS); /* s <- u^2+v^2 */
XCOPY(XU, XHP); /* h' <- u */
XCOPY(XV, XPP); /* p' <- v */
}
}
else
{
XSQR(Modulus, XS);
XCOPY(XS, MSqd);
}
if (HowMany == 1)
{
CalcMerit(XS, 2, Modulus, XMerit);
}
else
{
CalcMerit(XS, 2, MSqd, XMerit);
}
PrintSprectalData(XS, XMerit, 2);
if (MaxDimen > 2)
{
short I, J, K, NxtDimen;
USHORT *Umat, *Vmat;
USHORT R11[NE], R12[NE], R21[NE], R22[NE];
USHORT Vij[NE], Vjj[NE];
USHORT XA[NE], XB[NE];
USHORT *X, *Y, *Z;
/* ----------------- */
/* Initialize Arrays */
/* ----------------- */
INIT_VEC(X, MaxDimen);
INIT_VEC(Y, MaxDimen);
INIT_VEC(Z, MaxDimen);
INIT_MAT(Umat, MaxDimen);
INIT_MAT(Vmat, MaxDimen);
/* -------------------------------- */
/* Initialize U and V Matrices Etc. */
/* -------------------------------- */
if (HowMany == 1)
{
INIT_UMAT(Umat, MaxDimen);
INIT_VMAT(Vmat, MaxDimen);
}
else
{
XCOPY(Mult, XA);
XCOPY(Mult + NE, XB);
INIT_MAT_2(Umat, Modulus, MaxDimen);
INIT_MAT_2(Vmat, EONE, MaxDimen);
INIT_R_MAT(R11, R12, R21, R22);
}
MATPRT(Umat, 2, MaxDimen, "Initial U Matrix");
MATPRT(Vmat, 2, MaxDimen, "Initial V Matrix");
/* --------------------- */
/* S4: Advance Dimension */
/* --------------------- */
for (NxtDimen = 3; NxtDimen <= MaxDimen; ++NxtDimen)
{
USHORT Utt[NE];
INIT_VNEW(Vmat, NxtDimen, Modulus, MaxDimen);
if (HowMany == 1)
{
XMULT(XR, Mult, Tmp1); /* r * a */
XMOD(Tmp1, Modulus, XR); /* r <- (ra)mod m */
INIT_UNEW(Umat, NxtDimen, XR, MaxDimen);
FINISH_UV(Umat, Vmat, NxtDimen, Modulus, MaxDimen, XR);
}
else
{
SET_NEXT_RMAT(R11, R12, R21, R22, XA, XB, Modulus);
SET_NEXT_UROW(Umat, NxtDimen, R12, R22, MaxDimen);
COMPLETE_UV(Umat, Vmat, NxtDimen, Modulus, MaxDimen,
R12, R22);
}
MATPRT(Umat, NxtDimen, MaxDimen, "U Matrix At Step S5");
MATPRT(Vmat, NxtDimen, MaxDimen, "V Matrix At Step S5");
DOT(NxtDimen, Umat, NxtDimen, Umat, NxtDimen, MaxDimen, Utt);
XMIN(Utt, XS, XS);
J = 1, K = NxtDimen;
do
{
/* ------------- */
/* S5: Transform */
/* ------------- */
DOT(NxtDimen, Vmat, J, Vmat, J, MaxDimen, Vjj);
for (I = 1; I <= NxtDimen; ++I)
{
if (I != J)
{
DOT(NxtDimen, Vmat, I, Vmat, J, MaxDimen, Vij);
XCOPY(Vij, Tmp1);
XABS(Tmp1);
XRMULT(Tmp1, ETWO);
if (XGT(Tmp1, Vjj))
{
TRANSFORM(Vmat, I, Umat, J, Vij, Vjj,
NxtDimen, MaxDimen);
K = J;
}
}
}
/* --------------------- */
/* S6: Examine New Bound */
/* --------------------- */
if (K == J)
{
USHORT Ujj[NE];
DOT(NxtDimen, Umat, J, Umat, J, MaxDimen, Ujj);
XMIN(XS, Ujj, XS);
}
/* ------------- */
/* S7: Advance J */
/* ------------- */
if (J == NxtDimen)
{
J = 1;
}
else
{
++J;
}
}
while (J != K);
MATPRT(Umat, NxtDimen, MaxDimen, "U Matrix At End of S7");
MATPRT(Vmat, NxtDimen, MaxDimen, "V Matrix At End of S7");
/* ----------------------- */
/* S8: Prepare for Search */
/* ----------------------- */
CLR_VECTOR(NxtDimen, X);
CLR_VECTOR(NxtDimen, Y);
for (J = 1; J <= NxtDimen; ++J)
{
XDIV(XS, Modulus, Tmp1);
DOT(NxtDimen, Vmat, J, Vmat, J, MaxDimen, Z + XVL(J));
XRDIV(Z + XVL(J), Modulus);
XRMULT(Z + XVL(J), Tmp1);
XFLOOR(Z + XVL(J), Z + XVL(J));
XSQRT(Z + XVL(J), Z + XVL(J));
XFLOOR(Z + XVL(J), Z + XVL(J));
}
VECPRT(Z, NxtDimen, "Z Vector at Step S8");
/* ------------------------------------------------ */
/* S9 - S11 - The for control statement is Step S11 */
/* ------------------------------------------------ */
for (K = NxtDimen; K >= 1; --K)
{
/* -------------- */
/* S9: Advance Xk */
/* -------------- */
if (!XEQUAL(X + XVL(K), Z + XVL(K)))
{
XRADD(X + XVL(K), EONE);
INCREASE_Y(Y, Umat, K, NxtDimen, MaxDimen);
VECPRT(Y, NxtDimen, "Increased Y: ");
do
{
/* -------------- */
/* S10: Advance K */
/* -------------- */
++K;
if (K > NxtDimen)
{
XdotProd(NxtDimen, Y, Y, Tmp1);
XMIN(XS, Tmp1, XS);
}
else
{
XCOPY(Z + XVL(K), X + XVL(K));
XNEG(X + XVL(K));
DECREASE_Y(Y, Umat, K, Z + XVL(K),
NxtDimen, MaxDimen);
VECPRT(Y, NxtDimen, "Decreased Y: ");
}
}
while (K <= NxtDimen); /* end of Step S10 */
} /* end of if - Step S9 */
} /* end of for - Steps S9-S11 */
if (HowMany == 1)
{
CalcMerit(XS, NxtDimen, Modulus, XMerit);
}
else
{
CalcMerit(XS, NxtDimen, MSqd, XMerit);
}
PrintSprectalData(XS, XMerit, NxtDimen);
} /* end of for - Step S4 */
}
}
}